Computer tomographic method and a computer tomograph

ABSTRACT

In conventional backprojection, each pixel value is calculated separately. With N pixels N 3  additions are necessary. According to the invention, the summation is instead divided into partial sums which are then used to form new partial sums in a hierarchial way. A substantial computational reduction is obtained because a partial sum has approximately the same value for many neighboring points, and that interpolations are feasible. Instead of N 3  additions, the number of additions can be reduced to the order of magnitude of N 2  log N.

INTRODUCTION

The present invention regards computer tomography and methods ofevaluating tomography data. The method of the invention can be appliedto all kinds of tomographic imaging, whether registered by X-ray orradioactive absorption, by magnetic resonance (MR) or positron emissiontomography (PET). In all those cases are obtained experimental values inthe form of line integrals along multitudes of straight lines. Thetypical CT measurement consists of mounting an X-ray tibe on a rotatablecircular ring and set of detectors in a row opposite the tube.

The invention therefore regards a computer tomography method of the kindrecited in the preamble of claim 1.

In computerized tomographic imaging (CT) the inverse problem consists ofdetermining a function from its line integrals. We examine threedifferent types of measured data sets, with different scanninggeometries, namely two versions of parallel scanning (sinogram,linogram) and fan-beam scanning (even cone-beam scanning). We will studythe Filtered Backprojection Method (FBM) which is one of the mostwell-known and widely used methods. The reason for the widh use of thefiltered backprojection algorithm is the high image quality.

A serious problem in prior art evaluation of tomographic data in orderto obtain an image is the heavy computational complexity for theback-projection operator, which is of the order N³. It is the generalobject of the invention to obtain an improved method of transformingtomographic data into images, which is less complex and can giveimprovement either in computational speed or in an increased number ofvoxels/pixels in an image.

This object and other objects and advantages are obtained, according tothe invention, by means of a computer tomography method of the kindrecited, which is characterized by the steps recited in claim 1.Advantageous embodiments are recited in dependent claims.

In conventional backprojection each pixel value is calculatedseparately, which means that O(N) operations are required for eachpixel. Then the total number of operations will be of the order O(N³).Instead of directly doing the N³ additions the basic idea behind the newalgorithm is to divide the summation into partial sums, which then areused to form new partial sums in a dyadic and hierarchial way, seeFIG. 1. The reason for this procedure is the benefit in computationalreduction due to the fact that a partial sum has approximately the samevalue for many neighbouring points. Thus each value of a partial sum iscalculated once, but used several times when forming a new partial sumfrom two partial sums in a preceding level of the hierarchy.

Although in the present description, the invention is exemplified withforming hierachial partial sums in a dyadic way, which is presentlybelieved to be the best mode, it is also possible to sum more than twogroups, namely, each time a number which is a prime number. It is evennot necessary to use the same prime number in each stage.

In a preferred embodiment, families of lines are combined two and two,indicating that the total number of line families is a power of two.With a number of pixels/voxels of 512×512 an an equal number offamilies, i.e. angles of registration, this condition is fulfilled.However, without using the full advantages of the invention, it ispossible in the partial sum calculation to combine three (or five, orany prime number) space distributions by making sums for those pointswhere the lines for the largest angular difference intersect, and addinginterpolation values for the further distributions.

It is also possible to add line integrals for intermediate lines, forwhich no measurement was made, i.e. by setting intermediate(non-measured) line integrals to the mean value of the adjacent lines,for which a measurement was made.

Although the method is most easily understood for a case where all theline families belong to the same plane, which is the most commonconfiguration in practice, and which is exemplified as the fan-type, thesinogram and the linogram, it may also be adapted to the cone-beam case,where radiation from a source fills a cone, an detection is made in atwo-dimensional arry of detectors. Each detector then corresponds to aline through the field swept, and the contraption and the object aremutually rotated. As shall be shown, a considerable saving incalculations may be obtained also in this case.

In a special case, where each family of lines consists of lines whichare mutually parallel, it is possible, in a very favourable way, toobtain values for points which are equidistant in a rectangular screenpattern. This is obtained by making the two last line families, forwhich sums at intersections are calculated, be a mutually right angles.Thereby, in principle, the last repeated second step is performed by twocombined groups, symmetrized in relation to a direction at right anglesto an angle belonging to a family of lines which is at an extreme end ofsaid enumerable order, thus obtaining a set of sum values whichcorrespond to pixels in a cartesian system of screen pattern.

LIST OF FIGURES

FIG. 1 The summation process for N=8 . . . 21

FIG. 2 Division of the (k*,l*)-plane for I*(0,i,k*,l*) . . . 21

FIG. 3 Division of the (k*,l*)-plane for I*(1,i,k*,l*) . . . 22

FIG. 4 Division of the (k*,l*)-plane for I*(j,i,k*,l*) . . . 23

FIG. 5 Cone-beam geometry . . . 24

FIG. 6 Points of intersection for the partial sum I(j,i,k,l,m) . . . 25

2 THE FILTERED BACKPROJECTION ALGORITHM FOR SINOGRAM

In the parallel scanning geometry a cross-section of the human body isscanned by a set of equally spaced parallel X-ray beams for a number ofequally distributed directions. Thus the recorded detector valuescorrespond to a discrete set of values from the Radon transform. The(n-dimensional) Radon transform of a function f on R^(n) is defined by##EQU1## which is the integral of f over the hyperplane whose normal isθ.di-elect cons.S^(n-1), the unit sphere in R^(n), and with (signed)distance sεR¹ from the origin. The term sinogram is chosen since anattenuation function consisting of a δ-function, Dirac function, gives aRadon transform supported on a sinusoid in parameter space, i.e.,(θ,s)-space.

Before giving the inversion formula we will describe the backprojection.The (n-dimensional) backprojection operator R^(#) maps a function g onS^(n-1) ×R¹ into the set of functions on R^(n). The backprojectionoperator is defined by

    (R.sup.# g)(X)=∫.sub.Sn-1 g(θ,X·θ)dθ, X .di-elect cons.R.sup.n.

R^(#) is the dual operator of R. Hence

    ∫.sub.Sn-1 ∫.sub.R1 (Rf)(θ,s)g(θ, s)dsdθ=∫.sub.Rn ƒ(X)(R.sup.# g)(X)dX.

The filtered backprojection algorithm is a numerical implementation ofthe inversion formula

    ƒ=1/2(2π).sup.1-n R.sup.#H.sup.n-1 (Rf).sup.(n-1),

where H denotes the Hilbert transform and (·).sup.(n-1) denotes the(n-1)st derivative both with respect to the variable s.

3 Fast backprojection for sinogram

In this section we will develope a method for backprojection, requiringonly O(N² log N) operations. This method can be used on the(n-dimensional) backprojection operator but in order to describe themain idea it suffices to study the case n=2. Then ##EQU2## where X=##EQU3## ε]R² and θ= ##EQU4## ε S¹, or, with slightly differentnotation, ##EQU5##

Strictly we should write g(θ(v), x cos v+y sin v). We assume that g(θ,s)satisfies the condition g(-θ,-s)=g(θ,s) which implies thatg(v+π,-s)=g(v,s). Then ##EQU6## Let g(v,s) be samples at the points(g(v_(m),s_(n)) where ##EQU7## N=2^(J), J positive integer. In order tocompute the backprojection, we need to calculate the sums ##EQU8## fork² +l² ≦ ##EQU9## where k and l are integers.

Since g is defined on a rectangular grid we might use linearinterpolation for the variable s in the sum of formula 3.1 We chooseinstead the following procedure.

First extend the definition of g by increasing the number of samplingpoints in the s-coordinate with a factor e, using linear interpolation,and introduce the notation ##EQU10## for m=0, . . . , N-1 and ##EQU11##Then, take ##EQU12## where [a] denotes a rounded to the nearest integer.

Now, instead of directly doing the N³ additions we will divide thesummation into partial sums, which then are used to form new partialsums in a hierarchial way. The reason for this procedure is the benefitin computational reduction, due to the fact that a partial sum has thesame value for many neighbouring points. First we define ##EQU13## for0≦j≦J, 1≦i≦2^(J-j). Then we have the following simple relation

    I(j,i,k,l)=I(j-1,2i-1,k,l)+I(j-2i,k,l).

In the following description we will show that the number of points(k,l) where the value of the partial sum I(j,i,k,l) has to be calculatedis approximately twice as large as the number of points where the valueof a partial sum I(j-1,i,k,l) has to be calculated. Hence, since thenumber of partial sums with different i-values at level j is N/2^(j),the total number of summation points for all partial sums at level j isapproximately the same as for level j-1 and since there are J=² log Nj-levels we conclude that the total number of operations for all levelswill be proportional to N² log N.

The variable i just stands for a rotation about the origin. Thus it willbe convenient to make the following definition ##EQU14## Then we have

    I(j,i,k,l)=I*(j,i,k cos v(.sub.i-1)2j +l sin v.sub.(i-1)2j, l cos v(.sub.i-1)2j -k sin v.sub.(i-1)2j)                       (3.2)

and

    I*(j,i,k*,l*)=I*(j-1,2i-1,k*,l*)+I*(j-1,2,i,k* cos v.sub.2j-1 +l* sin v.sub.2j-1, l* cos v.sub.2j-1 -k* sin v.sub.2j-1).        (3.3)

B definition we have

    I*(0,i,k*,l*)=G(i-1,[ek*]).

We now see that I*(0,i,k₁ *,l₁ *)=I*(0,i,k₂ *,l₂ *) if the points (k₁*,l₁ *) and (k₂ *,l₂ *) are located in the same strip, see FIG. 2.

Now consider the next level, j=1. We have

    I*(1,i,k*,l*)=G(2i-2,[ek*])+G(2i-1,[e(k* cos v.sub.1 +l* sin v.sub.1)]),

and we conclude that I*(1,i,k₁ *,l₁ *)=I*(1,i,k₂ *,l₂ *) if (k₁ *,l₁ *)and (k₂ *,l₂ *) are located in the same region, see FIG. 3. Let##EQU15##

According to FIG. 3 it will be necessary to calculate I*(1,i,k*,l*) onlyfor a and b integers. Therefore we introduce, in general k* and l* willnot be integers,

    I(1,i,a,b)=I*(1,i,k*,l*)

where ##EQU16## Hence

    I(1,i,a,b)=G(2i-2,a)+G(2i-1,b)

Can we do a similar division of the (k*,l*)-plane for I*(j,i,k*,l*),j>1? Let ##EQU17## If two points (k₁ *,l₁ *) and (k₂ *,l₂ *) are in thesame region, see FIG. 4, we want I*(j,i,k₁ *,l₁ *)=I*(j,i,k₂ *,l₂ *),i.e. ##EQU18## for m=0, . . . , 2^(j) -1. In FIG. 4 we see that allregions approximately satisfy conditions (3.4) if ##EQU19## whichcorresponds to j≦J-1. In the same way as for j=1 we introduce

    I(j,i,a,b)=I*(j,i,k*,l*)                                   (3.5)

where ##EQU20## By equations (3.3) and (3.5) we see that we have thefollowing relation

    I(j,i,a,b)≈I(j-1,2i-1,a,b.sub.1)+I(j-1,2i, a.sub.1,b)(3.6)

where ##EQU21## and ##EQU22## By symmetry we have

    I(j,i,b,a)≈I(j-1,2i-1,b,a.sub.1)+I(j-1,2i,b.sub.1,a) I(j,i,-a,-b)≈I(j-1,2i-1,-a,-b.sub.1)+I(j-1,2i,-a.sub.1,-b)

and

    I(j,a,a,a)≈I(j-1,2i-1,a,a.sub.1)+I(j-1,2i,a.sub.1,a)

where ##EQU23## which means that we only have to compute a₁ and b₁ for0≦a≦|b|.

The procedure suggested by this discussion of the partial sums leads tothe following algorithm. We start by adding "extra" sample values bylinear interpolation. Then we form partial sums of adjacent pairwiseprojections (A projection m consists of all sample points (G(m,n) wherem is fixed.). The next step will be to form partial sums of fourprojections by adding two adjacent sums from the previous step.Thereafter we will subsequently form partial sums of 8,16,32, . . .projections. The iteration stops when full backprojection is obtained,that is, when the sum consists of all projections. A partial sum is onlydetermined in points where two lines intersect whose angles are v₀ andv_(2j-1), with respect to the (k*,l*)-plane, and whose perpendiculardistances to the origin are ##EQU24## and ##EQU25## respectively, n₁ andn₂ integers. If e≧2 the process is stopped before we obtainI(J-1,i,a,b). If the total number of points where the partial sums, atlevel j, are calculated is of the same order as the number of pixels,the process is stopped and the pixels are interpolated directly from thevalues of the partial sums at level j-1. The pixel values are thencalculated from equation (3.2). When e=1 we calculate all the partialsums up to I(J-1,i,a,b). Although the number of points whereI(J-1,i,a,b) is calculated is of the same order as the number of pixelsthey do not coincide. This is caused by ##EQU26## Therefore we have tocalculate the pixel values by equation (3.2) which is unnecessary if weinstead do the following. We want to calculate ##EQU27## and ##EQU28##But now ##EQU29## so the points where we calculate the sums are in factthe pixels. One problem is that each sum, (3.9) and (3.10), consist of2^(J-1) +1 elements. However in order for the summation process to workthe number of elements in a sum has to be a power of two. Therefore wecalculate ##EQU30## and ##EQU31## and ##EQU32## Note that since thepartial sums ##EQU33## and ##EQU34## both contain the projection m=N/2and the projection m=0˜m=N, we do the following. We assign theprojections the values G(N/2,n).tbd.0 and G(0,n).tbd.0 in one of thepartial sums.

The total number of different values for ##EQU35## is ¢N√2.

If j and i are fixed then the number of points where we calculateI(j,i,a,b) is approximately the area of the circle, ##EQU36## divided bythe are of one region (rhomb), ##EQU37## Hence the total number ofpoints at level j, ##EQU38## will be ##EQU39## When j is small(j=1,2,3.), the overlapping of the circle by rhombs is not so efficientbut then ##EQU40## Hence the total number of points for all levels is##EQU41## for e≧2 and ##EQU42## for e=1. 4 The filtered backprojectionalgorithm for linogram

In 1987 a new reconstruction algorithm was developed by Edholm [Edh87].It is called the linogram method and is similar to the direct Fouriermethod but it is based on a different scanning geometry. In the parallelscanning geometry for sinogram, the projection lines, X-rays, are takenfor a number of equally distributed directions which requres thedetectors to rotate around the object. The linogram projections actuallyconsist of two sets of projections. In the first linogram the detectorsand the transmitters lie equally spaced on fixed lines parallel to thex-axis and the projection angles are -45°≦θ≦45°. Instead of equallydistributed directions, δθ=constant, we now get δtanθ=constant. In thesecond linogram the detectors and transmitters lie equally spaced onfixed lines parallel to the y-axis and 45°≦θ≦135°, (δcot θ=constant).

Instead of the Radon transform it now will be natural to study the twolinogram operators R₁ and R₂. ##EQU43## The term linogram is chosensince an attenuation function consisting of a δ-function gives a R₁ -and R₂ -transform supported on a line in parameter space, i.e.,(u,v)-space.

There exist dual operators, backprojection operators, R₂ ^(#) and R₂^(#) given by ##EQU44## We have ##EQU45## The filtered backprojectionalgorithm for linogram is a numerical implementation of the inversionformula ##EQU46## where H denotes the Hilbert transform and (·).sup.(i)denotes the derivative both with respect to the variable u.

5 Fast backprojection for linogram

Here the procedure of the summation for the backprojection is analogousto the one used for sinogram. It suffices to study the backprojectionR^(#) ₁ and projection data from R₁ since R^(#) ₂ and R₂ corresponds toa rotation about the origin of 90° degrees. We want to determine (R^(#)₁ g)(x,y) in the pixels (x,y)= ##EQU47## N=2^(J) where J is a positiveinteger and k, l integers belonging to the disk k² +l² ≦ ##EQU48## Thenthe function g(u,v) has to be sampled in the points g(u_(m),v_(n)),##EQU49## m integer, for all g(u_(m),v_(n)) where the correspondingprojection line intersect the unit circle. Therefore the number ofparallel projection lines depends on the projection angle. For instanceif n=0 we have |m| ≦ ##EQU50## but if n= ##EQU51## we have |m| ≦##EQU52## We will divide the summation. ##EQU53## into partial sums inthe same way as we did for sinogram. In order for the summation processto work, the number of elements in a sum has to be a power of two.Therefore we start by splitting the summation in the following way##EQU54## Then we divide the sum ##EQU55## into partial sums (Theprocedure for ##EQU56## is similar). Defin ##EQU57## for 0≦j≦J-1,1≦i≦2^(J-j-1). The procedure suggested by this discussion of the partialsums is analogous to the method described in section 3 on parge 8. Apartial sum I(j,i,k,l) is only determined in points where two projectionlines intersect, one projection line corresponds to v.sub.(v-1).2j+1 andthe other to v_(i).2j. If the projection line, which corresponds to(u_(a),v.sub.(i-1).2j+1), intersects the projection line, whichcorresponds to (u_(b),v_(i).2j), at the point (k_(j),i,a,b,l_(j),i,a,b)in the (k,l)-plane we will use the notation

    I(j,i,a,b)=I(j,i,k.sub.j,i,a,b,l.sub.j,i,a,b).

The projection line x=u_(m) -v_(n).y is the straight line through thepoints (m+n, -N/2) and (m-n, N/2) in the (k,l)-plane. We have ##EQU58##Since we only are interested in points on the disc ##EQU59## we get

    |b-a|≦2.sup.j -1.

We obtain a similar relation for linogram as relation (3.6) forsinogram. We have

    I(j,i,a,b)≈I(j-2i-1,a,b.sub.1)+I(j-2i,a.sub.1,b)

where ##EQU60## and ##EQU61## One major advantage of linogram comparedto sinogram is that in the equations for a₁ and b₁ the terms which notare integers depend only on the difference b-a. Then if we use linearinterpolation the interpolation factors will be constant ifb=b-a=constant. Hence if we introduce the notation

    I(j,i,a,b)=I(j,i,a,a+b)

we finally obtain the simple and very effective summation algorithm

    I(j,i,a,b)≈(1-c)I(j-1,2i-1,a,b-A)+cI×(j-1,2i-1,a,b-A-1)

    +(1-c)I(j-1,2i,a+A,b-A)+cI(j-1,2i,a+A+1,b-A-1)

where ##EQU62## =A+c. A integer and 0≦c<1. Since c=0 for b=0,±(2^(j) -1)and by using that we have the same interpolation factors for ±b, we onlyhave to compute ##EQU63## for 1≦b≦2^(j) -2 (Let 1≦b≦2^(j) -2 and##EQU64## =Ac, A integer and 0<c<1. Then ##EQU65## =-A-1+(1-c). ). Forsinogram, however, we had to compute a₁ and b₁ for all combinations of(a,b), a≧b. Therefore the total number of times a₁ and b₁ have to becompute, at level j, is approximately 2(2^(j) -1)N for sinogram. Hencealthough the number of summations points (a,b) is approximately the samefor linogram and sinogram the total computational complexity for a₁ andb₁ is O(N) for linogram but O(N²) for sinogram.

Finally we calculation R^(#) ₁ g according to equation (5.11). Let

    I1(a,b)=I(J-1,1,a,b)

and I2(a,b) for the corresponding partial sum ##EQU66## where acorresponds to the direction v₋₁ and b=a+b corresponds to the directionv_(-N/2). Then ##EQU67## for 0≦l≦N/2.

The total computational complexity for the filtered fast backprojectionalgorithm for the linogram geometry is the same as for the sinogramgeometry, i.e., N² log N.

6 Cone-beam and Fan-beam geometry

The most common scanning geometries in medical tomography are thecone-beam geometry in the three-dimensional case, and the analogousfan-beam geometry in the two-dimensional case. In cone-beam scanning anX-ray source is rotating around the object along a circle while emittingcones of X-rays, which are recorded by a two-dimensional detector. Infan-beam scanning the cone-beams are screened off to a single plane andconsequently recorded by a one-dimensional detector. Feldkamp, Davis andKreiss [Fel84] have developed an approximate 3D-reconstruction algorithmfor cone-beam projections. The computational complexity is of the orderO(N⁴) for an image ƒ of N×N×N× pixels. By using the same technique as wedid for the previous backprojection algorithms we will reduce the numberof operations to O(N³ log N).

Next we will describe the cone-beam geometry. Let the position of theX-ray source be denoted by S(λ), ##EQU68## where 0≦λ<2π and r constant.Hence the X-ray source rotates around the object along a circle in thexy-plane. The detector plane is always perpendicular to the vector S(λ),see FIG. 5. Hence it will be practical to introduce the plane spanned bythe normed vectors p(λ) and z, where p(λ) is the vector in the xy-planewhich is perpendicular to S(λ), see FIG. 5. Then the measured data canbe represented by a discrete set of values from ##EQU69## where L(λ,p,q)is the straight line through the points S(λ) and D(p,q)=pp(λ)+qz. Nowthe approximate inversion formula can be written as [Fel84] ##EQU70##where ##EQU71## and ##EQU72## i.e., we multiply the projection datag(λ,p,q) with a weight function ##EQU73## and then perform arampfiltering in the variable p. The final backprojection in (6.12) iscomputationally the most heavy part of the inversion algorithm. Themajor difference between this backprojection and all the previous is theweight function A(λ,x). Now let g(λ,p,q) be sampled at the pointsg(λ_(a),p_(b),q_(c)) where r>1 and N=2^(J), with J integer, ##EQU74##and we want to estimate ƒ(x) in the points ##EQU75## k,l and m integers,for ##EQU76## Thus ##EQU77## where ##EQU78## and

    x.sub.k,l,m ·p(λ.sub.a)=x.sub.k sinλ.sub.a -y.sub.1 cosλ.sub.a =(x·p)(a,k,l).

If we now divide the sum into partial sums in the same way as we did insection 3 we get ##EQU79## Now if the interval is small, i.e. j small,we can make the approximation

    A(a,k,l)≈A(a.sub.j,i,k,l)                          (6.13)

where a_(j),i =i.2^(j) -(2^(j) +1)/2, the average angle. Thus ##EQU80##Define ##EQU81## Then as before

    I(j,i,k,l,m)=I(j-1,2i-1,k,l,m)+I(j-1,2i,k,l,m)

A partial sum is only determined in points where two projection lines,projected on the xy-plane, intersect one projection line coming fromS(λ.sub.(i-1)-2j) and the other from S(λ_(i).2j-1), see Fig 6. Now sinceA(λ,x) obviously can't be approximated by a constant for all λ, exceptx=0, the iteration process is stopped after a few steps, at level j₀,and the pixels are interpolated directly from the values of the partialsums at level j₀. We have ##EQU82##

At a first glance the method seems to be far less efficient for thisbackprojection than for the previous ones. But if we consider that thelevel j=1 reduces the number of summations by half and j=2 toapproximately a fourth of the original number, and so on, we see that itis almost as efficient as the other. The accuracy of the approximation(6.13) depends strongly on the distance between the point (k,l) and theorigin, i.e., √k² +l² . But since the inversion formula (6.12) iscorrect only for points in the plane z=0 and the artifacts in the otherplanes z=m increase with the absolute value of m the region of interestusually is √k² +l² +m² <<r. For (k,l) fixed point the accuracy of (6.13)depends on the angle difference λ_(i).2j-1 -λ.sub.(i-1).2j not on thenumber of elements, 2^(j), in the partial sum. ##EQU83##

Hence if N=2^(J) corresponds to that the iteration process is stopped atlevel j₀ then 2N corresponds to that the iteration process is stopped atlevel j₀ +1 in order to have the same accuracy for the approximation(6.13). Then if the region of interest consists of n³ pixels the numberof operations for the backprojection will be of the order O(n² N log N).In the two-dimensional case, fan-beam scanning, the procedure isidentical, m=0. Hence if the region of interest is of size n² thecomputational complexity will be O(nN log N).

References

[Edh87] P. R. Edholm. G. T. Herman, Linograms in image reconstructionfrom projections IEEE Transactions on medical imaging, Vol. MI-6, pp.301-307, 1987.

[Fel84] L. A. Feldkamp, L. C. Davis, J. W. Kress, Practical cone beamalgorithms Journal of optical soc. Am., Vol. A6, pp. 612-619, 1984.

[Nat86] F. Natterer, The mathematics of computerized tomography JohnWiley & Sons Ltd, Chichester, 1986.

[Rad17] Uber die bestimmung von funktionen durch ihre integralwerte lagsgewisser mannigfaltigkeiten Ber. Verh. Sachs. Akad. 69, p. 262, 1917.

Instead of using the approximation (6.13) we can even perform themultiplication of the weight function in an iterative way. In order todescribe the main idea is suffices to study the two-dimensional case,z≡0, the extension to the case z≠0 is obvious. We introduce thenotations x= ##EQU84## εR², L(λ,p)=L(λ,p,0) and ##EQU85## Letx(j,i,p₁,p₂) denote the point of intersection between the linesL(λ_(i-1))·2j, 2p₁ /N) and L(λ_(i)·2j-1, 2p₂ /N). Then we will use thenotation

    I(j,i,p.sub.1,p.sub.2)=I(j,i,x(j,i,p.sub.1,p.sub.2)).

We obtain a similar relation as (3.6)

    I(j,i,p.sub.1,p.sub.2)≈I(j-1,2i-1,p.sub.1,d.sub.1)+I(j-1,2i, d.sub.2,p.sub.2),

where d₁ and d₂ are calculated in a similar way as equations (3.7) and(3.8) for a₁ and b₁. In general d₁ and d₂ will not be integers. In orderto obtain the correct value of the weight function we will apply thefollowing procedure. We will only consider the partial sumI(j-1,2i-1,p₁,d₁) the procedure for I(j-1,2i,d₂,p₂) is analogous. Westart by only consider the term which corresponds to a=(i-1).2^(j)+2^(j-1) -1, i.e., ##EQU86## where x₁ =x(j-1,2i-1,p₁,d₁), x₅₇₂ 2=x(j-1,2i-1,p₁,[d₁ ]), x₃ =x(j-1,2i-1,p₁,[d₁ ]+1), d₁ =[d₁ ]+c₂, 0≦c₂ <1and c₁ +c₂ =1. Therefore we will use the following approximation##EQU87## where a=(i-1)·2^(j) +2^(j-1) -1. In order for theapproximation to be accurate for all terms in the partial sumI(j-1,2i-1,p₁,d₁) all weight functions which belong to the partial sumhave to satisfy the condition ##EQU88## for (k-1).2^(j)≦a1,a2≦(i-1).2^(j) +2^(j-1) -1 and [d₁ ]≦d≦[d₁ ]+1. The error obtainsits maximum when a1=(i-1).2^(j) and a2=(i-1).2^(j) +2^(j-1) -1. Itsuffices to consider the case ##EQU89## for i=1 and λ=2^(j-1) -1. Thepoints x(j-1, 1,p₁,d₁) and x(j-1,1,p₁,d) are located on the line L(0,2p₁/N) which can be written ##EQU90## Let x(j-1,1,p₁,d₁)=x(t₁) and(j-1,1,p₁,d)=x(t). Then ##EQU91## We have ##EQU92##

The number of points of intersection increase with λ. In fact we have|t₁ -t|∝ ##EQU93## Hence ##EQU94## for some constant c. At the firststeps of the iteration the error will be extremely small ∝1/N² and thenincrease to ∝1/N at the final step.

Now we return to the three-dimensional case, z≠0. We introduce thenotations ##EQU95## and

    I(j,i,p.sub.1,p.sub.2,z)=I(j,i,x(j,i,p.sub.1,p.sub.2),z).

We have the relation

    I(j,i,p.sub.1,p.sub.2,z)≈I(j-1,2i-1,p.sub.1,d.sub.1 z)+I(j-1,2i,d.sub.2,p.sub.2,z),

We will use the same approximation as before the only difference is thatwe have to adjust the value of the variable z. ##EQU96## wherea=(i-1)·2^(j) +2^(j-1) -1, λ=(i-1)·2^(j) +(2^(j-1) -1)/2 the averageangle and

    z·A(λ,x.sub.1)=z.sub.1 ·A(λ,x.sub.2)=z.sub.2 ·A(λ,x.sub.3).

In general z₁ and z₂ are not integers but since the weight function andthe points of intersection x₂ and x₃ are independent of z we can performan interpolation of the variable z in the usual way. The accuracy of theapproximatation in the z-direction depends on the previous condition,i.e., ##EQU97## for (i-1)·2^(j) ≦a≦(i-1)·2^(j) +2^(j-1) -1.

What is claimed is:
 1. A computer tomography method for obtainingamapping of a variable in an object, said variable being measured in theform of line integral values over a multitude of lines, crossing saidobject, the multitude of measured lines being subdivided into familiesof lines, which families are enumerable in an order determined bysuccessive relative angular shifts between patterns of lines insuccessive families, characterized in thata) in a first step, groups arecombined out of said families, each group consisting of a first primenumber of families of lines having adjacent position in said order, saidlines of said highest and lowest order of families forming a multitudeof intersection points between lines belonging to one famly and linesbelonging to another family, and calculating for each of saidintersection points of the sum of the line integral values for the tworespective said lines intersecting at that point, and further addingvalues interpolated to said points from eventual intermediate families,said first step resulting, for each said group of families, in aplurality of points in a plane, for which points individual sum valuesare stored, said sums forming space distributions, one for each saidgroup, b) in a successive step combining a further prime number ofpreviously combined groups, with their calculated space distributions,representative of adjacent positions in said order into further groupsrepresentative of a number of families which is said further primenumber times that of an immediately preceding step, having successivepositions in said order, and calculating for points of intersectionbetween those lines belonging to the families of the lowest and highestorder within the number of families comprised in each combined group, ofsums of values representative of the two respective space distributions,and summing thereto values for the aforesaid points of intersectiontaken for distributions belonging to eventual intermediate groups,either for points coinciding with those having calculated individualsums using such sums, or, for points not so coinciding, adding valuesinterpolated from adjacent points in a respective space distribution,forming in said successive step higher order space distributions, to anumber which is reduced by the factor of said further number compared tothat of the starting point of said successive step, c) and repeatingsaid second step until obtension of a final space distribution.
 2. Acomputer tomography method of claim 1, characterized in that saidinterpolations are performed by linear interpolation of values belongingto points of intersection regarding a line participating in a respectiveintersection point, for which an interpolated value is desired.
 3. Acomputer tomography method of claim 1, characterized in that saidinterpolations are performed taking into account neighbour distributionpoints in the plane, belonging to more than one line.
 4. A computertomography method of claim 1, characterized in that said first primenumber is
 2. 5. A computer tomography method of claim 4, characterizedin that each said further prime number is
 2. 6. A computer tomographymethod of claim 5, wherein each family of lines consists of lines whichare mutually parallel, characterized in that the last repeated secondstep is performed by two combined groups, symmetrized in relation to adirection at right angles to an angle belonging to a family of lineswhich is at an extreme end of said enumerable order, thus obtaining aset of sum values which correspond to pixels in a cartesian system ofscreen pattern.
 7. A computer tomography having a plurality of detectorsfor sensing intensities along a multitude of lines, some of which areintersecting each other, said lines intersecting an object to betomographed, said plurality of detectors being connected to a calculatorfor calculating pixel values for the intersected object, characterizedin that the said calculator is provided with means for calculatingpartial sums according to claim 1.